      SUBROUTINE GFLUXI(NLL,TLEV,NW,DW,DTAU,TAUCUM,W0,COSBAR,UBARI,    &
                        RSF,BTOP,BSURF,FTOPUP,FMIDP,FMIDM)

!     Legacy Mars GCM v24
!     Mars Climate Modeling Center
!     NASA Ames Research Center

!  THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITIONS
!  FOR THE INFRARED FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
!  THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
!  MEASURED FROM THE TOP OF EACH LAYER.  THE TOP OF EACH LAYER HAS  
!  OPTICAL DEPTH ZERO.  IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
!  HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
!  FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
! THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY 
! VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER
!
! NLL            = NUMBER OF LEVELS (NLAYERS + 1) MUST BE LESS THAT NL (101)
! TLEV(L_LEVELS) = ARRAY OF TEMPERATURES AT GCM LEVELS
! WAVEN          = WAVELENGTH FOR THE COMPUTATION
! DW             = WAVENUMBER INTERVAL
! DTAU(NLAYER)   = ARRAY OPTICAL DEPTH OF THE LAYERS
! W0(NLEVEL)     = SINGLE SCATTERING ALBEDO
! COSBAR(NLEVEL) = ASYMMETRY FACTORS, 0=ISOTROPIC
! UBARI          = AVERAGE ANGLE, MUST BE EQUAL TO 0.5 IN IR
! RSF            = SURFACE REFLECTANCE
! BTOP           = UPPER BOUNDARY CONDITION ON IR INTENSITY (NOT FLUX)
! BSURF          = SURFACE EMISSION = (1-RSFI)*PLANCK, INTENSITY (NOT FLUX)
! FP(NLEVEL)     = UPWARD FLUX AT LEVELS
! FM(NLEVEL)     = DOWNWARD FLUX AT LEVELS
! FMIDP(NLAYER)  = UPWARD FLUX AT LAYER MIDPOINTS
! FMIDM(NLAYER)  = DOWNWARD FLUX AT LAYER MIDPOINTS
!
!----------------------------------------------------------------------C

      use grid_h
      use constants_h, only: PI
      use radinc_h
      use radcommon_h, only: planckir

      implicit none

      integer, parameter :: NLP = 101    ! MUST BE LARGER THAN NLEVEL 

      INTEGER :: NLL, NLAYER, L, NW, NT, NT2
      REAL*8  :: TERM, CPMID, CMMID
      REAL*8  :: PLANCK
      REAL*8  :: EM,EP
      REAL*8  :: COSBAR(L_NLAYRAD), W0(L_NLAYRAD), DTAU(L_NLAYRAD)
      REAL*8  :: TAUCUM(L_LEVELS), DTAUK
      REAL*8  :: TLEV(L_LEVELS)
      REAL*8  :: WAVEN, DW, UBARI, RSF
      REAL*8  :: BTOP, BSURF, FMIDP(L_NLAYRAD), FMIDM(L_NLAYRAD)
      REAL*8  :: B0(NLP),B1(NLP),ALPHA(NLP),LAMDA(NLP),XK1(NLP),XK2(NLP)
      REAL*8  :: GAMA(NLP),CP(NLP),CM(NLP),CPM1(NLP),CMM1(NLP),E1(NLP)
      REAL*8  :: E2(NLP),E3(NLP),E4(NLP)

      REAL*8  :: FTOPUP, FLUXUP, FLUXDN
      real*8  :: TAUMAX = L_TAUMAX 

!======================================================================C
 
!     WE GO WITH THE HEMISPHERIC CONSTANT APPROACH IN THE INFRARED
 
      IF (NLL .GT. NLP) STOP 'PARAMETER NL TOO SMALL IN GLUFV'

      NLAYER = L_NLAYRAD

      DO L=1,L_NLAYRAD-1
        ALPHA(L) = SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) )
        LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI

        NT2   = TLEV(2*L+2)*10.0D0-499
        NT    = TLEV(2*L)*10.0D0-499

        B1(L) = (PLANCKIR(NW,NT2)-PLANCKIR(NW,NT))/DTAU(L)
        B0(L) = PLANCKIR(NW,NT)
      END DO

!     Take care of special lower layer

      L        = L_NLAYRAD
      ALPHA(L) = SQRT( (1.0-W0(L))/(1.0-W0(L)*COSBAR(L)) )
      LAMDA(L) = ALPHA(L)*(1.0-W0(L)*COSBAR(L))/UBARI

      NT    = TLEV(2*L+1)*10.0D0-499
      NT2   = TLEV(2*L)*10.0D0-499
      B1(L) = (PLANCKIR(NW,NT)-PLANCKIR(NW,NT2))/DTAU(L)
      B0(L) = PLANCKIR(NW,NT2)

      DO L=1,L_NLAYRAD
        GAMA(L) = (1.0-ALPHA(L))/(1.0+ALPHA(L))
        TERM    = UBARI/(1.0-W0(L)*COSBAR(L))

!       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
!       BOTTOM OF THE LAYER.  THAT IS AT DTAU OPTICAL DEPTH

        CP(L) = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM
        CM(L) = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM

!       CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
!       AT THE TOP OF THE LAYER, THAT IS ZERO OPTICAL DEPTH

        CPM1(L) = B0(L)+B1(L)*TERM
        CMM1(L) = B0(L)-B1(L)*TERM
      END DO
 
!     NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
!     FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
!     WARNING IF DTAU(J) IS GREATER THAN ABOUT 35 (VAX)
!     WE CLIP IT TO AVOID OVERFLOW.

      DO L=1,L_NLAYRAD

!       CLIP THE EXPONENTIAL HERE.

        EP    = EXP( MIN((LAMDA(L)*DTAU(L)),TAUMAX))
        EM    = 1.0/EP
        E1(L) = EP+GAMA(L)*EM
        E2(L) = EP-GAMA(L)*EM
        E3(L) = GAMA(L)*EP+EM
        E4(L) = GAMA(L)*EP-EM
      END DO
 
!     B81=BTOP  ! RENAME BEFORE CALLING DSOLVER - used to be to set
!     B82=BSURF ! them to real*8 - but now everything is real*8
!     R81=RSF   ! so this may not be necessary

!     Double precision tridiagonal solver

      CALL DSOLVER(NLAYER,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP,       &
                   BSURF,RSF,XK1,XK2)
 
!     NOW WE CALCULATE THE FLUXES AT THE MIDPOINTS OF THE LAYERS.
 
      DO L=1,L_NLAYRAD-1
        DTAUK = TAUCUM(2*L+1)-TAUCUM(2*L)
        EP    = EXP(MIN(LAMDA(L)*DTAUK,TAUMAX)) ! CLIPPED EXPONENTIAL 
        EM    = 1.0/EP
        TERM  = UBARI/(1.-W0(L)*COSBAR(L))

!       CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
!       BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH

        CPMID    = B0(L)+B1(L)*DTAUK +B1(L)*TERM
        CMMID    = B0(L)+B1(L)*DTAUK -B1(L)*TERM
        FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
        FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID

!       FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT

        FMIDP(L) = FMIDP(L)*PI
        FMIDM(L) = FMIDM(L)*PI
      END DO
 
!     And now, for the special bottom layer

      L    = L_NLAYRAD

      EP   = EXP(MIN((LAMDA(L)*DTAU(L)),TAUMAX)) ! CLIPPED EXPONENTIAL 
      EM   = 1.0/EP
      TERM = UBARI/(1.-W0(L)*COSBAR(L))

!     CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
!     BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH

      CPMID    = B0(L)+B1(L)*DTAU(L) +B1(L)*TERM
      CMMID    = B0(L)+B1(L)*DTAU(L) -B1(L)*TERM
      FMIDP(L) = XK1(L)*EP + GAMA(L)*XK2(L)*EM + CPMID
      FMIDM(L) = XK1(L)*EP*GAMA(L) + XK2(L)*EM + CMMID
 
!     FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT

      FMIDP(L) = FMIDP(L)*PI
      FMIDM(L) = FMIDM(L)*PI

!     FLUX AT THE TOP LEVEL

      EP   = 1.0
      EM   = 1.0
      TERM = UBARI/(1.0-W0(1)*COSBAR(1))

!     CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
!     BOTTOM OF THE LAYER.  THAT IS AT DTAU  OPTICAL DEPTH

      CPMID  = B0(1)+B1(1)*TERM
      CMMID  = B0(1)-B1(1)*TERM

      FLUXUP = XK1(1)*EP + GAMA(1)*XK2(1)*EM + CPMID
      FLUXDN = XK1(1)*EP*GAMA(1) + XK2(1)*EM + CMMID
 
!     FOR FLUX WE INTEGRATE OVER THE HEMISPHERE TREATING INTENSITY CONSTANT

      FTOPUP = (FLUXUP-FLUXDN)*PI
 
      RETURN
      end subroutine gfluxi
